High— order jamming crossovers and density anomalies 
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We demonstrate the existence of high-order jamming crossovers in systems of particles with re- 
pulsive contact interactions, which originate from the collapse of successive coordination shells. At 
zero temperature, these crossovers induce an anomalous behavior of the bulk modulus, which varies 
non-monotonically with the density, while at finite temperature they induce density anomalies con- 
sisting in an increased diffusivity upon isothermal compression and in a negative thermal expansion 
coefficient. We rationalize the dependence of these crossovers on the softness of the interaction po- 
tential, and relate the jamming crossovers and the anomalous diffusivity through the investigation 
of the vibrational spectrum. 
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Repulsive potentials with a finite range are investi- 
gated in a variety of different contexts, including the 
jamming transition of non thermal systems, and density 
anomalies of complex liquids. In the first case, these 
potentials are considered as models for the interaction 
between soft particle systems such as emulsions, bub- 
ble rafts and granular materials, which are not affected 
by Brownian motion. These systems undergo a jam- 
ming transition marking the onset of mechanical rigidity 
when the volume fraction crosses a threshold <fij. See 
Refs. [H, HJ for recent reviews. At the jamming transition 
each particle is forced to touch particles in its first coor- 
dination shell, and the mean contact number per particle 
jumps to the isostatic value Z- lso , which is the minimum 
value required for mechanical stability. Above 4>j the ex- 
cess contact number AZ = Z — Z- lso grows as a power 
law in j , and is related to a length scale diverging at 
the transition as AZ -1 , known as isostatic length. The 
jamming transition is characterized by an abundance of 
soft vibrational modes: the density of states D(tu) satis- 
fies Debye scaling up to a characteristic frequency scaling 
as ui* oc AZ, above which D(ui) flattens. 

Similar potentials, and in particular repulsive poten- 
tials with a cutoff distance D and with a smaller hard- 
core radius, are also commonly investigated at finite tem- 
perature to reproduce the density anomalies of water and 
of some other network-forming liquids. See Ref. Q for a 
recent review. A density anomaly is a region of the phase 
diagram where the diffusivity increases upon isothermal 
compression, which usually overlaps with a region char- 
acterized by a negative thermal expansion coefficient. In 
network-forming liquids, these anomalies are related to 
soft vibrational modes characterized by a rigid rotation 
of tetrahedral structures, known as rigid unit modes JiJ. 
Even though an abundance of soft modes characterizes 
both the jamming transition and the density anomalies, 
these two phenomena have not been related before, possi- 
bly because density anomalies occur at volume fractions 
well above </>j [5j-|9[. 



In this Letter we show that the jamming transition 
is the first of a series of high-order jamming crossovers. 
These occur on increasing the volume fraction as particles 
are forced to make contact with those of subsequent co- 
ordination shells. The geometrical signatures of the jam- 
ming crossovers are oscillations in the rate of formation 
of new contacts on compression. The mechanical ones 
include an anomalous volume fraction dependence of the 
elasticity of the system, whereby the bulk and shear mod- 
uli vary non-monotonically with volume fraction. We 
show that density anomalies are the finite temperature 
counterpart of the jamming crossovers, and clarify the 
relation between these phenomena via the study of the 
soft vibrational modes. 

Model - We perform Molecular Dynamics simulations 
of 50:50 binary mixtures of particles with diameter Di = 
1 and D s = Di/IA and mass M = 1. Two particles 
i and j interact if they have a positive overlap, 5ij = 
Dij — Tij > 0, where Dij is their average diameter, and 
Tij their distance. When this is the case the interaction 
potential is 



V(n 3 ) 



h (ta- 
ct \Di 



(1) 



where the parameter a controls the softness of the in- 
teraction potential, larger values of a corresponding to 
softer potentials. M, Di and k = 1 arc our units 
of mass, length and energy. Finite temperature sim- 
ulations have been performed in the NVT and NPT 
ensembles, considering systems of N = 10 3 particles. 
Zero temperature properties have been investigated by 
quenching random configurations of N = 10 4 parti- 
cles to the closest energy minimum via the conjugate- 
gradient algorithm [l(| ■ We have considered a large vol- 
ume fraction range, varying from the jamming thresh- 
old </> ~ <t>j (<j)j ~ 0.84 in 2d, <j)j ~ 0.64 in 3d) up 
to = 3, and considered the following values of a: 
a = 1.25,1.5 (Hertzian), 1.75, 2 (Harmonic), 2.5 and 3. 
The a = 1 case corresponds to the Jagla potential with 
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FIG. 1: (color online) (a) Volume fraction dependence of the 
diffusion coefficient averaged over the two species, for different 
interaction potentials. For each potential, the temperature is 
chosen so that the minimum value of D in the range of <j) 
studied here is D ~ fCP 6 . (b) Thermal expansion coefficient. 



no hard-core repulsion [5j|. We report in the following 
results of 2d systems; analogous results hold in 3d. 

T > properties - We start by showing that sys- 
tems of particles interacting via purely repulsive finite 
range potentials are not characterized by a single den- 
sity anomaly as previously reported , but by a 
series of successive density anomalies. Indeed, as illus- 
trated in Fig. [IJi, there exist successive volume fraction 
ranges in which the diffusivity increases upon isothermal 
compression. An analogous result has been observed in 
models of polymer stars Likewise, Fig. [TJj shows 
the existence of multiple volume fraction ranges charac- 
terized by a negative thermal expansion coefficient. The 
diffusivity anomalies depend on the interaction poten- 
tial: while for soft potentials (large a) a single anomaly 
is observed, more anomalies are observed on increasing 
the stiffness of the interaction. Likewise, the potential 
influences the number and the strength of the anomalies 
in the thermal expansion coefficient. In the following, 
we clarify the microscopic origin of these anomalies, and 
show that they are the finite temperature counterparts 
of high-order jamming crossovers. 

T = properties - We have investigated the vol- 
ume fraction dependence of geometrical and mechanical 
properties of jammed configurations in a large volume 
fraction range. To this end we have considered volume 
fractions with a spacing of A<f> = 5 • ICP 2 far from the 
jamming transition and A<f) = 1 CP 2 close to it; data are 
averaged over 50 independent jammed configurations for 
each value of </>. Measurements taken at volume fractions 
differing by A<f> have been used to estimate volume frac- 
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FIG. 2: (color online) Volume fraction dependence of the 
mean contact number Z (a), of its volume fraction derivative 
dZ/d<f> (b), of the bulk modulus K (c), of the non affinity pa- 
rameter x (d), and of the normalized characteristic frequency 
zu (e). Each column refers to a different a value, as indicated. 



tion derivatives numerically [13j. Figs.[5k,b illustrate the 
volume fraction dependence of the mean contact num- 
ber and its first derivative. At the jamming transition, 
Z jumps to the isostatic value, while at higher volume 
fractions it grows with superimposed oscillations. Con- 
sistently, dZ/d<j) diverges at the jamming transition, and 
then varies in an oscillating manner. As the divergence of 
dZ/dcf) at the jamming transition corresponds to the col- 
lapse of the first coordination shell, the successive volume 
fraction ranges in which dZ/d(f> increases can be inter- 
preted as high-order jamming crossovers, during which 
each particle gradually makes contact with neighbors in 
successive coordination shells. In line with the fact that 
the shell structure is lost at high density, the crossovers 
become smoother as the volume fraction increases, and 
eventually Z increases linearly with (f>. As an aside we 
note that for soft potentials there is a volume fraction 
range in which the average contact number is constant, 
Z = 6; here all particles in contact are Voronoi neighbors, 
and the value Z = 6 is fixed by Euler's theorem for pla- 
nar graphs. The jamming crossovers also influence me- 
chanical properties such as the bulk modulus K, whose 
volume fraction dependence is show Fig. The effect 
of the crossovers on K depend on a. Indeed, minima and 
maxima of K and dZ/d(j> are closely correlated for small 
a, and become less correlated as a increases, indicating 
that the smaller a, the greater the influence of the new 
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FIG. 3: (color online) Panels a,b,c: overlap probability dis- 
tribution function Pz{$) for different values of the volume 
fraction separated by A(j> — 5 x 10 -2 , as indicated. Data 
are shifted vertically for clarity. A negative 5 corresponds to 
the distance between the surfaces of non-interacting parti- 
cles. The inset of panel (d) shows that the intersection vol- 
ume fraction 4> n is to a good approximation independent on 
the softness of the interaction a, specified in the legend. 



contacts on the properties of the system. We finally note 
that the results of Fig. [5] do not support the conjecture 
of a transition at a volume fraction 4>d > 4>J 14 1- 



We elucidate the microscopic evolution of the system 
at the jamming crossovers investigating the volume frac- 
tion dependence of the overlap p.d.f. Pz($), normalized 
across the interacting particles so that f Q Pz(S)d5 = Z. 
For monodisperse particles Pz is fully equivalent to the 
radial distribution function. Since dZ/d<f> increases when 
particles start interacting with those of a new shell, we 
show in Fig. [3] Pz (§) at volume fractions immediately 
following that of the first minimum of dZ/d(j>, for dif- 
ferent values of a. For a = 3 the making of new con- 
tacts does not strongly affect the distribution. Pz{&) is 
simply shifted towards larger values of 6, which suggests 
the presence of a mostly afBne deformation. For a = 2, 
on the other hand, the formation of new contacts leads 
to a transformation of Pz(8), whose net outcome is the 
creation of a new peak at small positive <5. This peak 
shifts towards positive S values on further compressing 
the system. At smaller a, Pz($) has a diverging peak 
at 5 — > + , and the transformation of the distribution 
driven by contact formation leads to an increase of this 
peak. The higher order jamming crossovers lead to anal- 
ogous changes in Pz{5). These results clarify that the 
jamming crossovers lead to a rearrangement of the force 
network of the system, which is more pronounced the 
smaller a. In particular, we note that that on increasing 
the volume fraction the shell structure observed in Pz (5) 
for 6 < 0, is transformed into a shell structure for pos- 
itive overlaps for large a, while conversely it is washed 



out for smaller a. 

To explain the role of the softness of the interaction 
potential in the shape of the overlap probability distri- 
bution, we introduce the concept of intersection volume 
fraction <p n . We define </> n as the sum of the volumes 
of intersection of any pair of interacting particles, nor- 
malized by the total volume of the system. Close to the 
jamming transition, where the intersection volumes do 
not overlap, <j> = <pj + <pm so <Pn is linearly dependent 
on <p and potential independent. In addition, here <pn{<P) 
scales as the q-th moment of Pz(8), (f>n oc (5 q ), as the 
volume of intersection of two particles grows as 6 q for 
small 5, where q = 1.5 for d = 2 (q = 2 for d = 3). 
Fig. [3H (inset) shows that 4> n is also approximately po- 
tential independent at higher volume fractions, so we can 
regard is as being fixed by <j> throughout. The energy of 
the system, on the other hand, scales as the ev-th mo- 
ment of Pz{8), E oc (<5 a ). Accordingly, constant-;/) en- 
ergy minimization protocols search for a Pz distribution 
that minimizes (5 a ) subject to the constraint of constant 
(S q ), or to related constraints away from jamming. The 
effect of a on the outcome can be understood by consid- 
ering the energy e(<5,) = (<5f + S")/a of two overlaps Si 
and Sj related by the constraint of constant intersection 
volume fraction, 8f + 5 1 = const. Then for a > q it 
is energetically more favourable for the system to make 
the overlaps equal, while for a < q the energy is mini- 
mized by making them maximally different. This finding 
explains the abundance of small contacts in the overlap 
distribution for small a, and the qualitative change in Pz 
as a crosses q. 

The a dependence of the overlap distribution indi- 
cates that, the stiffer the interaction potential, the more 
heterogeneous the structure of the system. This sug- 
gests an increase of the non-affine response of the sys- 
tem on decreasing a. We estimate the degree of affinity 
of the system by comparing the actual bulk modulus K 
with that computed in the affine (Born) approximation, 
-/tag. If the system responds affinely to the compres- 
sion, then K = K a ff, otherwise K < K a g. Accordingly, 
the strength of the non-affine response, i.e. the relevance 
of the fluctuation term of the stress tensor [l5[, can be 
quantified via the parameter 



K. 
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ail 



K 



0<X< 1. 



Then x — when the response is affine, while x ~ * 1 
when the response is highly non-affine. Fig. confirms 
the expectation that stiffer potentials give rise to a less 
affine response, as the typical value of x decreases on in- 
creasing q. In addition, the comparison of Figs. andHJl 
reveals the presence of a clear anticorrelation between 
X and K , which indicates that the jamming crossovers 
induce an increase of the degree of non-affinity of the 
system, much like the jamming transition. We also note 
that the oscillations of xi'P) suggest that the non-affine 
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correlation length that diverges at the jamming transi- 
tion [HQ is a non-monotonic function of the volume 
fraction. 

Connecting T = and T > The mechanical 
anomalies observed at T = 0, and the dynamical ones 
observed at T > 0, have a closely related volume frac- 
tion dependence. This suggests that they have the same 
physical origin. We elucidate this mechanism by exploit- 
ing recent results correlating the thermal [n| and shear 
17, [IK induced relaxation of glassy systems, and their 
vibrational modes. These studies have shown that the 
relaxation proceeds through localized events, which oc- 
cur where the low frequency modes are quasi-localized. 
In this picture, low-frequency modes are defects allowing 
for the relaxation of disordered particle systems (T^]. In 
particular, since the activation energy of a soft mode is 
correlated with its eigenfrequency [19j, one expects the 
dynamics to speed up when the typical frequency of the 
soft modes decreases. We have checked this expectation 
by studying the volume fraction dependence of the av- 
erage eigenfrequency w of the lowest modes (5%), nor- 
malized by the average eigenfrequency of all modes. Re- 
sults are illustrated in Fig. for different interaction 
potentials, w exhibits oscillations as the volume frac- 
tion increases, which are anti-correlated with those of 
the diffusion coefficient. In addition, w is also strongly 
anti-correlated with the parameter \- These results con- 
firm the important role of the soft modes in determining 
both the relaxation dynamics as well as the mechanical 
properties of the system. 

Conclusions - We have demonstrated the existence 
of high-order jamming crossovers in soft particle sys- 
tems with increasing density. These occur when parti- 
cles start to touch neighbors in subsequent coordination 
shells. The mechanical manifestation of these crossovers 
is the anomalous behavior of elastic properties such as 
the bulk and the shear (data not snown) moduli, which 
vary non-monotonically with the volume fraction. Their 
dynamical manifestation are density anomalies, which in- 
clude both an increased diffusivity on compression and 
a negative thermal expansion coefficient. We have re- 
lated the T = and the T > anomalies to the emer- 
gence of soft vibrational modes in correspondence with 
the crossovers, confir ming their important role in the re- 
laxation dynamics 1||- 18 1 . These results suggest that the 



scaling relation between relaxation time and vibrational 
dynamics found in normal liquids po| may also hold for 
anomalous liquids. We finally note that the anomalies are 
observed at the high-order jamming crossovers, and not 
at the jamming transition, where the diffusion coefficient 
decreases monotonically on compression. One may spec- 
ulate that this is so because at the jamming crossovers 
there is a coexistence of compressed bonds, and of new 
almost uncompressed bonds. These bonds may play the 
role of the two interaction length scales which are known 
to induce dynamical anomalies in models of water Q, 



and are responsible for the rigid unit modes leading to a 
negative thermal expansion coefficient in network form- 
ing systems These two length scales are not present 
at the jamming transition. In this respect, it would be 
interesting to investigate the spatial structure of the soft 
modes at the jamming crossover. 
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